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Abstract 

An algorithm for computing {2, 3}, {2, 4}, {1, 2, 3}, {1, 2, 4}-inverses 
and the Moore-Penrose inverse of a given rational matrix A is established. 
Classes A{2, 3} s and yl{2,4} s are characterized in terms of matrix prod- 
ucts (R*A) 1 R* and T*(AT*)\ where R and T are rational matrices with 
appropriate dimensions and corresponding rank. The proposed algorithm 
is based on these general representations and the Cholesky factorization 
of symmetric positive matrices. The algorithm is implemented in pro- 
gramming languages MATHEMATICA and DELPHI and illustrated via ex- 
amples. Numerical results of the algorithm, corresponding to the Moore- 
Penrose inverse, are compared with corresponding results obtained by 
several known methods for computing the Moore-Penrose inverse. 

AMS Subj. Class.: 15A09, 68Q40. 
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1 Introduction 

Let C be the set of complex numbers, C mx ™ be the set of m x n complex 
matrices, and C™ x ™ is a subset of C mxn consisting matrices of rank r: C™ x ™ = 
{X 6 c mx ™ | rank(A) = r}. As usual, C(x) denotes the set of rational functions 
with complex coefficients in the variable x. The set of m x n matrices with 
elements belonging to C(x) is denoted by C(x) mxn . By I r and / we denote 
the identity matrix of the order r, and identity matrix of an appropriate order, 
respectively. By O is denoted an appropriate null matrix. 

For any matrix A of the order mxn consider the following matrix equations 
in X, where * denotes conjugate and transpose: 

(1) AXA = A (2) XAX = X (3) (AX)* = AX (4) (XA)*=XA. 



'Corresponding author 
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In the case m = n we also consider equations 

(5) AX = XA (l fc ) A k+1 X = A k . 

For a sequence S of elements from the set {1, 2, 3, 4, 5, l fc }, the set of matrices 
obeying the equations with corresponding indicative numbers contained in iS is 
denoted by ^4{<S}- A matrix from j4{5} is called an S- inverse of A. The matrix 
X = A^ is said to be the Moore-Penrose inverse of A satisfies equations (l)-(4). 
The group inverse A# is the unique {1,2,5} inverse of A, and exists if and only 
if ind(A) = min{fc : rank(A fc+1 ) = rank(A fe )} = 1. A matrix X = A D is said 

to be the Drazin inverse of A if (l fe ) (for some positive integer k), (2) and (5) 
are satisfied. In the case ind(^4) = 1, the Drazin inverse of A is equal to the 
group inverse of A. If A is nonsingular, it is easily seen that ind(A) = and 
A D =A-\ 

The rank of generalized inverse X is important, and it will be convenient to 
consider the subset A{i,j, k} s of A{i,j, k}, consisting {i, j, fc}-inverses of rank 
s (see [1]). 

In the literature are known various methods for computing the Moorc- 
Pcnrose inverse (see for example [1], [24]). The most commonly implemented 
method in programming languages is the Singular Value Decomposition (SVD) 
method, that is implemented, for example, in the "pinv" function from Matlab, 
as well as in the standard MATHEMATICA function " Pseudolnverse" [26]. This 
method is very accurate, but time consuming when the matrix is large. Other 
well-known methods are Greville's algorithm, the full rank QR factorization by 
Gram-Schmidt orthonormalization (GSO), and iterative methods of various or- 
ders [1]. A number of expansions of the Moore-Penrose inverse can also be used 
to develop direct methods [15], [21]. 

A class of direct methods for computing pseudoinverses is derived from the 
full-rank factorization A = PQ of m x n matrix A of rank r, where P is m x r, 
Q is r x n, and P, Q are both of rank r. These methods are investigated in 
many papers (see for example [1, 16, 21, 24]). After the full-rank factorization, 
we have the general representation of the Moore-Penrose inverse A^ = P\ 
where Q t = Q*{QQ*)- 1 ,P^ = (P*P)- 1 P*. General representations for various 
classes of {2}-inverses and the Drazin inverse are obtained in [21]. 

Chen et all derived a deterministic iterative algorithm for computing the 
Moore-Penrose inverse and rank of matrix A e C mx ™ in [4]. This algorithm 
is called successive matrix powering and it is based on successive squaring of 

" P Q " 



a composite matrix T = 



O I 



where P = (I - [3 A* A), Q = (3 A* and (3 



is a relaxation parameter. Wei established successive squaring algorithm to 
approximate the Drazin inverse in [25]. The Drazin inverse is expressed in 

r P 

the form of successive squaring of the composite matrix T= 



P=(I - l3A k+1 ),Q = (3A k . 



O I 



where 
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In the paper [5] , Courrieu proposed an algorithm for fast computation of the 
Moore-Penrose inverse of real matrices, which is based on known reverse order 
law (eq. 3.2 from [15]), and on the full rank Cholesky factorization of possibly 
singular, symmetric positive matrices (Theorem 4 from [6]). 

In the present paper we use the L [/-factorization from [6]. An arbitrary 
matrix A has an LU-factorization if it can be expressed as the product A = LU 
of a lower-triangular matrix L and an upper triangular matrix U. When it 
is possible, we say that A has an LU- decomposition. It turns out that this 
factorization (when exists) is not unique. If L has l's on it's main diagonal, 
then it is called a Doolittle factorization. If U has l's on its main diagonal, 
then it is called a Crout factorization. When L = U* , it is called the Cholesky 
decomposition. In each of these cases, the following is valid: 

a* = u*tf = u*{uu*)- 1 {l*l)- 1 l*. 

An implementation of the Cholesky factorizations in MATHEMATICA can be 
found on the web site 

http : I /math. fullerton.edu/mathews/n2003/CholeskyBib. html. 

This paper is a generalization of the paper [5] to sets of {2, 3}, {2, 4}-inverses 
and to the set of rational matrices. 

Many numerical algorithms for computing the Moore-Penrose inverse lack 
numerical stability. Also, when rounding error is present, we have to identify 
some small quantity as being zero. Moreover, it is well-known that the Moore- 
Penrose inverse is not necessarily a continuous function of the elements of the 
matrix. The existence of this discontinuity is an additional problem in the 
pseudoinverse computation. It is clear that cumulative round off errors should 
be totally eliminated. This is possible only by symbolic computation. During 
the symbolic implementation, variables are stored in the "exact" form or can 
be left "unassigned" (without numerical values), resulting in no loss of accuracy 
during the calculation [10]. 

Algorithms for computing generalized inverses of polynomial and/or rational 
matrices are so far based upon the Leverrier-Faddeev algorithm and the Grevile's 
algorithm. Computation of the Moore-Penrose inverse of polynomial and/or 
rational matrices which uses the Leverrier-Faddeev algorithm is investigated in 
[7, 9, 10, 11, 23]. An algorithm of the Leverrier-Faddeev type for computing 
the Moore-Penrose inverse of a polynomial matrix is introduced in the paper 
[10]. Implementation of this algorithm, in the symbolic computational language 
MAPLE, is described in [9]. Furthermore, in [9] it is described an implementation 
of the algorithm for computing the Moore-Penrose inverse of a singular rational 
matrix. 

A representation and corresponding algorithm for computing the Drazin 
inverse of a singular one-variable polynomial matrix of arbitrary degree are 
introduced in [8], [18]. Corresponding algorithm for two- variable polynomial 
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matrix and its implementation is introduced in [2]. Also, an effective version of 
given algorithm is established in the paper [2] . 

A general finite algorithm for computing various classes of generalized in- 
verses of a polynomial matrix is introduced in [20]. This algorithm is based on 
the Leverrier-Faddeev algorithm. 

Computation of the Moore- Penrose inverse of one- variable polynomial and/or 
rational matrices, arising from the Grevile's algorithm, is introduced in [17]. 
Corresponding two-dimensional case is investigated in [22]. 

The Moore- Penrose inverse is used in the evaluation of the least square solu- 
tion of linear system Ax = b, even with rank deficient matrices [1]. In fact, the 
Moore-Penrose inverse is defined as that matrix which, when postmultiplicd 
by b, yields the minimum-length least-square solution x of the possibly incon- 
sistent equation Ax « 6, for any b. Also, the Moore-Penrose inverse can have 
valuable applications in neurocomputational learning procedures [5]. Moreover, 
in the literature it is known a number of applications of generalized inverses of 
polynomial matrices [9, 10, 11, 12, 13, 14]. 

This paper is a first attempt to compute k} generalized inverses of 
one- variable rational matrices using the method from [5] . 

In the second section we characterize classes A{2,3} s , A{2,4} s , A{1,2,3} 
and A{1,2,4} in terms of matrix products (R*A^R* and T*(AT*)\ where 
R and T are rational matrices with appropriate dimensions and correspond- 
ing rank. Using these representations, we introduce a method for computing 
fc}-inverses of prescribed rank s of a given rational matrix A. When A is 
a constant matrix, in two partial cases (R = A or T = A), we get an algorithm 
for computing the Moore-Penrose inverse, alternative with corresponding one 
introduced in [5]. 

Algorithm introduced in this paper is implemented in programming package 
MATHEMATICA, and it is applicable to rational and constant matrices. Corre- 
sponding algorithm, applicable only to constant matrices, is also implemented 
in the programming language DELPHI. Symbolic implementation in MATHEMAT- 
ICA is illustrated via examples in Section 3. We especially consider the partial 
case of the implementation, which computes the Moore-Penrose inverse of a con- 
stant matrix. This partial case of the implementation is compared with several 
known methods for computing the Moore-Penrose inverse. 

2 Representations of {i,j,k} inverses for rational 
matrices 

In the following lemma we modify known representations for {2,3}, {2,4}- 
inverses of prescribed rank, introduced in [1]. We also extend these representa- 
tions, known for complex matrices, to the set of one-variable rational matrices. 

Lemma 2.1 Let A E C(x)" ixn and < s < r, mi,ni > s be chosen integers. 
Then the following general representations for pseudoinverses are valid: 
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(a) A{2,4} s = {(YAYY\ Y G C(x)" lXm , Y A G C(x)" lX ™}. 

(b) A{2,3} s = {Z(AZ)1\ Z G C(2;)" xmi , G C^)™ 1 } . 

Proo/. (a) The inclusion A{2,4} S D{ (YA^Y \ YeC(x) n ^ xm , YAeC(x) r ^ xn } 
can be proved in a similar way as in [1]. 

To prove the opposite inclusion, choose an arbitrary X G A{2, 4} s . Consider 
a full-rank factorization X = FG, FeC(i)f s ,GeC(i)f m . Since X G A{2}, 
we get 

FGAFG = FG 

or 

F(GAF - I S )G = O. 

This implies 

GAP = J a . 

Now, it is not difficult to verify F G (GA){1, 2, 3, 4}, or equivalently F = (GA)^ . 
Consequently, 

X = (GA)*G G {(YA) f Y\ Y G C(x) sxm , YA G C(x)f "} . 

Using 

{(yA)ty|y g c(i) sxm ,yiec(i)f n } c { (ya^y\y& c(x) ni xm , fag c(x)^ 1 xn } 

we prove part (a). 

Part (b) can be verified in a similar way. □ 

Remark 2.1 In the case m\=n\ = s, in the case of constant matrices, we get 
an improvement in the proof of Theorem 6 and Theorem 7 from [1] (p. 63). 

Analogous representations of {1, 2, 3} and {1, 2, 4}-inverses we derive in the 
case s = r = rank(A). 

Lemma 2.2 Let A G C(x)™ xrl and mi,ni > r be chosen integers. Then the 
following statements are valid for the sets A{1, 2, 4}, A{1, 2, 3} and the Moore- 
Penrose inverse: 

(a) A{1,2,4} = {(YA)^Y\ Y G C(a;)" lXm , YA G C(:r); 11 x ™} . 

(b) A{1,2,3} = {Z(AZY\ Z G C(a:) nxm \ G C(x)™ 1 } . 

(c) At = = A*{AA*y. 

Now we are in a position to propose the next theorem for computing {2, 3}, 
{2,4} inverses of prescribed rank as well as {1,2,3} and {1,2,4} inverses of a 
given matrix A G C(x)™ xn . This theorem is a customization of Lemma 2.1 to 
generalized LU factorization from [5] and [6] . 
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Theorem 2.1 Consider rational matrix A G C(x)" ixn . Let < s < r be ran- 
domly chosen integer and assume that m\,ni are positive integers satisfying 
mi,ni>s. Then the following statements are valid: 

(a) 

A{2,4} S = {L(L*L)- 2 L*(R*A)*R*\ Re C(x)™ xm , R*Ae C(x)^ 1 xn } , (2.1) 

where (R*A)*(R*A) = LL* is the Cholesky factorization and L* is without the 
zero rows. 

(b) 

A{2,3} S = {T*(AT*)*L(L*L)- 2 L*\ Te C(x)™ 1 x ™, AT* e C(x)™ 1 } , (2.2) 

where (AT*)(AT*)* = LL* is the Cholesky factorization and L* is without the 
zero rows. 

(c) 

A{1,2,4} = {L(L*L)- 2 L*(R*A)*R*\ ReC(x)^ lxn \ R* AeC{x)^ xn } , (2.3) 

where (R*A)*(R*A) = LL* is the Cholesky factorization and L* is without the 
zero rows. 

(d) 

A{1,2,3} = {T*(AT*)*L(L*L)- 2 L*\ TeC(x)" 11 xn ,AT* eC(x) 1 l nxmi } , (2.4) 

where (AT*)(AT*)* = LL* is the Cholesky factorization and L* is without the 
zero rows. 

(e) 

A^=L(L*L)- 2 L*(A*A)*A*, (2.5) 

where (A* A)* (A* A) = LL* is the Cholesky factorization and L* is without the 
zero rows, or 

A** = A*{AA*)*L{L*L)- 2 L*, (2.6) 

where (AA*)(AA*)* = LL* is the Cholesky factorization and L* is without the 
zero rows. 

Proof, (a) Various expressions for computing the Moore-Penrose inverse of the 
matrix product {AB)^ are considered in [15]. We use the following: 

(AB)t = B*(A* ABB*)^ A* . (2.7) 

Applying (2.7) in the case A = R*A, B = I, the Moore-Penrose inverse 
(R*A)^ can be found as 

(iTA) f = ((R*A)*(R*A))\R*A)*. (2.8) 

There is an unique upper triangular matrix S with exactly n — s zero rows, 
such that S*S — (R*A)*(R*A), where the computation of S is an application 
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of the extension of the usual Cholesky factorization from [5], [6] on matrix 
(R*A)*(R*A). Removing the zero rows from S, one obtains anrxn matrix of 
rank r, denoted by L*. The following is evident: 

(R*A)*(R*A) = S*S = LL*. (2.9) 

Applying (2.9) in (2.8), we get 

(i2M)t = {LL*)\R*Af. (2.10) 
Applying now (2.7) in the case A = L, B = L* , one can verify the following 

(LL*)t =L(L*L)- 1 {L*L)- 1 L*. (2.11) 

Multiplying (R*A)^ by R* from the right, in view of (2.10) and (2.11) we obtain 

(R*A)*R* = L(L*L)- 2 L*(R*A)*R*. 

Now, the proof follows from Lemma 2.1, part (a). 

(b) This part of theorem can be proved in a similar way as part (a) , applying 
part (b) from Lemma 2.1 and A = I, B = AT*. Also, in this case m and mi 
appears instead of n and m, respectively. 

Parts (c), (d) and (e) can be proved applying Lemma 2.2. □ 

Using Theorem 2.1, we now state the following algorithm which generates 
classes A{2,4} s and A{2,3} s . 

Algorithm 2.1 Choose mxn rational matrix A and consider randomly chosen 
mi x m rational matrix R, where m\ = m and m is arbitrary integer > r, or 
m = n and mi is arbitrary integer > r. 

Step 1. If n = ni then compute G :— (AR*)(AR*)* and set n — m and logical 
variable trans = True; 
else compute G := (R*A)*(R*A). 

Step 2. Find Cholesky factorization of matrix G — LL* and drop zero rows from 
L*. 

Step 3. If trans then return R*(AR*)*L(L*L)- 2 L* ; 
else return L(L* L)^ 2 L* (R* A)* R* . 

This algorithm is applicable to class of rational matrices if we implement 
them in symbolic programming languages like MATHEMATICA, MAPLE etc. Our 
implementation is developed in MATHEMATICA. However, because of the prob- 
lems with the simplification in rational expressions, this algorithm is not conve- 
nient for the implementation in high level programming languages such as C++, 
DELPHI, VISUAL BASIC etc. Therefore, our implementation in language DELPHI is 
applicable only for constant matrices. 
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3 Examples 

Example 3.1 In this example we consider constant matrices. Let A gCf 
and -ReC®* 6 be the following matrices: 



A = 
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) 




V o 


-1 








-2 


1 / 



Applying the function ModGinvCholesky [A,R] , described in Appendix, we ob- 
tain 






and 



A{l,2,A} 



( -ife 



V 



1 

10 

1 

13 




i 

2 
-1 
1 


1 

10 



3 
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Let us mention that conditions of Theorem 2.1, part (c) are valid. 

Example 3.2 Let us consider matrix A of rank 3: 
A={{x+l,x,5},{x+2,x,3},{x-l,x,l},{x+3,x,2}} . 

Choose the following matrix R of rank 2: 

R={{x+1 , 2} , {x+1 , 2} , {x+1 , 3} , {x+1 , 3}} . 

In accordance with part (a) of Theorem 2.1, function ModGinvCholesky [A, R] 

generates the following {2, 4} -inverse of A of rank 2: 
In[3] := ModGinvCholesky[A, R] 

Ontt'l] — ff -21-30 x+A x 2 -21-30 x+Ax 2 56+80 x~Ax 2 56+80 x-4x 2 i 

KJaL V°\ — 11 49+140 z+204 x 2 ' 49+140 z+204 x' 1 ' 49+140 z+204 x' 2 ' 49+140 z+204 x' 1 J' 



( -2x_(rf+2x) -2x (17+2 x) 2 x (43+2 x) 2 x (43+2 x) -| 

I 49+140 x+204 x' 2 ' 49+140 z+204 x' 1 ' 49+140 z+204 x 2 ' 49+140 z+204 x 2 J' 
r 14+34 z+40 3; 2 14+34 x+AO x 2 21+44 x+AO x 2 21+AAx+AQx 1 

I 49+140 a;+204 i 2 ' 49+140 a;+204 x 2 ' 49+140 a;+204a: ;4 ' 49+140 1+204 i 



^}} 



Example 3.3 In this example we choose matrices A andT satisfying conditions 
imposed in part (d) of Theorem 2.1. Then an {1,2, 3} -inverse is generated in 
the output: 

In[4] := A = {{1 + x,x,5, 2 + x,x,3}, {-1 + x,x, 1,3 + x,a;,2}, 

{— 2 + a;, a;, 1, 3 + x,x,2}, {— 3 + x, — 1 +a;, 1, 1 + x,x, 1}}; 

7n[5] := T = {{1 + x, 2, 2 + x, 1, -1 + x, 3}, {2 + cc, 3, 3 + x, 1, -2 + cc, 2}, 
{3 + a;, 3, 3 + a;, -1, -2 + a;, 1}, {2 + x, 3, 3 + x, 4, -1 + a;, 1}, 
{2 + x, 3, 3 + x, -1,-1 + x, 1}}; 

7n[6] := ModGinvCholesky[A,T] 
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r)lif\(i\— UCl 1 —1 ()\ f 5031-13465 a;+14101 ai 2 -130 a; 3 -975a: 4 
WUL[U\ ±, 1, UJ-, \ 30186^78744 a;+63024 x A +49340 a; 3 + 7800 a; 4 ' 

-70434+89855 a:-4908 a: 2 + 18453 a; 3 +6615 x 4 -75465+82542 a: + 12803 a: 2 +22101 a; 3 +6180 x l 
30186-78744 a;+63024 a; 2 +49340 a; a +7800 a; 4 ' 30186-78744 1+63024 aj 2 +49340 a; 3 + 7800 a; 4 ' 
-25155+20168 a;+11555 a; 2 +4153 a; 3 +540 x 4 f ( 5031-11455 a;+7726 a; 2 +5905 a; 3 +975 x 4 
30186-78744 1+63024 x'- 1 +49340 a; 3 + 7800 a; 4 J ' I 30186-78744 u+63024 a^+49340 a; 3 + 7800 a; 4 ' 
-70434+90587 a;+28674 a; 2 -3784 a; 3 -1695 x 4 75465-85296 a=-45272 a: 2 -4707 a; 3 +540 x 4 
30186-78744 a; + 63024 a; 2 +49340 a; a + 7800 a; 4 ' 30186-78744 aj + 63024 a: 2 +49340 a; 3 + 7800 a; 4 ' 
-25155 + 22250 x+19052 a: 2 +4011 a: 3 + 180 x 4 -| r 5031 - 15463 a= + 15407 a: 2 +8530 a: 3 +975 x 4 
30186-78744 a; + 63024 a; 2 +49340 a; 3 + 7800 a; 4 J ' I 30186-78744 a; + 63024 a; 2 +49340 a; 3 + 7800 x 4 ' 
-50310+38335 ai + 81884 ai 2 + 21697 a: 3 +735 a= 4 75465-86214 a:-56095 a: 2 + 1091 a: 3 +2780 x 4 
30186-78744 ai+63024 a: 2 +49340 a; 3 + 7800 a; 4 ' 30186-78744 a;+63024 aj 2 +49340 a; 3 +7800 a: 4 ' 
-35217+49192 ai+543 a: 2 -12483 a; 3 -2540 x 4 1 ( -12+3071 a=-13389 ai 2 +4760 a: 3 + 1950 x 4 
30186-78744 a;+63024 a; 2 +49340 a; 3 + 7800 a; 4 > ' I 30186-78744 a;+63024 a; 2 +49340 a; 3 + 7800 a; 4 ' 
-96960+246044 x+ 1 19747 a; 2 -42630 a; 3 - 15150 x 4 106974-261537 ai-153182 ai 2 + 21230 ai 3 + 11200 x 4 

30186-78744 a;+63024 a^+49340 a; 3 +7800 a; 4 ' 30186-78744 a;+63024 a^+49340 a; a + 7800 a: 4 ' 
-29838+61289 a:+78394 a: 2 +22290 a: 3 +2000 x 4 -j r 5031-17461 ai + 16713 a; 2 + 17190 a: 3 +2925 x 4 
30186-78744 a;+63024 a^+49340 a; a + 7800 a: 4 J ' I 30186-78744 a;+63024 a^+49340 a; a +7800 a; 4 ' 

-140868+87781 a=+221884 x 2 + 111187 a: 3 + 15885 x 4 -166023+97482 a;+251041 a; 2 + 118599 a; 3 + 16220 x 4 
30186-78744 1+63024^^+49340 a; a + 7800 a; 4 ' 30186-78744 1+63024:^+49340 a: a + 7800 a: 4 ' 

-65403+39808 ai+75665 a; 2 +28527 ai 3 +3260 x 4 1 1 
30186-78744 a;+63024 a; 2 +49340 a; a +7800 x 4 ' ' 

Example 3.4 In this example we generate {1,2, 4} -inverse using the following 
matrices A and R: 

A={{x+l ) x ) 5} ) {x+2 ) x ) 3} ) {x-l,x ) l},{x+3 ) x ) 2},{x-2,x,l} ) {x+3,x,2}}. 

R={{l+x , 2 , 2+x , 1 , -1+x} , {2+x , 3 , 3+x , 1 , -2+x} , {3+x , 3 , 3+x , -1 , -2+x} , 
{2+x , 3 , 3+x , 4 , -1+x} , {2+x , 3 , 3+x , - 1 , -1+x} , { 1+x , 2 , 2+x , 1 , -1+x}} . 

Jn[9] := ModGinvCholesky[A, R] 

n,,f\n]— ff 1596-2292 ai + 2542 a; 2 -5593 + 2166 a:+4766 x 2 2310-5520 a: + 1282 a: 2 

CMI[aj— \\ -41601 + 39942 a;-27634 a: 2 ' 83202-79884 rr + 55268 a: 2 ' -41601 + 39942 ai-27634 a; 2 ' 
13573-13626 a:+7944 a: 2 10549-4878 a=+7922 a: 2 1596-2292 z+2542 a: 2 -j 

41601-39942 1 + 27634^^ ' -83202 + 79884 a;-55268 x'- 1 ' -41601 + 39942 a;-27634 a;^ J ' 

f 23961-36414 x+38983 a: 2 -5084 a; 3 -38171+3108 a:+76654 a: 2 -9532 x 3 60564-14028 a;+34011 a: 2 + 2564ai 3 
I -83202 a; + 79884 a; 2 -55268 a; 3 ' 4a; (41601-39942 a: + 27634a: 2 ) ' 83202 a:-79884 a; 2 +55268 a; 3 ' 
74774-105294 a:+62993 a: 2 -15888 a: 3 29743-69888 a:-4194 a: 2 + 15844 a; 3 23961-36414 a:+38983 a: 2 -5084 a: 3 - 
83202 ai-79884 a; 2 +55268a; 3 ' 166404 a;-159768 a; 2 +110536 a; 3 ' -83202 aj + 79884 a; 2 -55268 a: 3 - 

f 16905-21399 a:+20869 a: 2 -7 (-5257+1862 a;+4414 a: 2 ) 18228+3465 ai+14261 a; 2 
I 83202-79884 a;+55268 a; 2 ' 4 (41601-39942 a;+27634 a: 2 ) ' -83202+79884 a;-55268a; 2 ' 

-7 (5446-5735 a;+2597 a; 2 ) 8281+25270 a:+12302 a: 2 16905-21399 a;+20869 a; 2 -| -| 
83202-79884 a;+55268 a; 2 ' 166404-159768 z + 110536 a; 2 ' 83202-79884 rr+55268 a; 2 J J 



Example 3.5 In this example we choose matrices A and T satisfying condi- 
tions imposed in part (b) of Theorem 2.1. Then an {2, 3}-inverse of rank 2 is 
generated: 



/n[10] 
/n[ll] 
7n[12] 



= A = {{x + l,x,5},{x + 2,x,3},{x - l,a;,l}, {x + 3, a;, 2}}; 
— T — {{x + 1, 2, x — 1}, {x + 2, l,x - 1}} 
= ModGinvCholesky[A,T\ 



r)a,f\T9] — ff 41-139 a; + 88 a; 2 +25 a; 3 +a; 4 30-99 a: + 83a; 2 +31 a= 3 +3a; 4 

UUL[LZ\ — \\ 329-1168 a; + 984a; 2 +380 a; 3 +35 a: 4 ' 329-1168 a; + 984 a; 2 +380 a; 3 +35 a; 4 ' 

55-239 a: + 222 a; 2 +97a: 3 +9 a; 4 85-290 a; + 228 a: 2 +82 a; 3 + 7 a: 4 -| r - 136 + 69 x + 207 x 2 +35 a; 3 +a; 4 
329-1168 a; + 984a: 2 +380a: 3 +35 a; 4 ' 329-1168 a; + 984 a; 2 + 380 a; 3 + 35 x 4 J ' 1 329-1168 a: + 984 a; 2 +380 a; 3 +35 a: 4 ' 
69-170a; + 40a: 2 +26a; 3 +3a: 4 -38+3 a: + 373 a: 2 + 117 a; 3 +9 a: 4 31 -254 a: + 246 a; 2 +82 a: 3 + 7 a; 4 -| 

329-1168 a; + 984a; 2 +380 a; 3 +35 a; 4 ' 329-1168 a: + 984 a; 2 + 380 a; 3 + 35 a: 4 ' 329-1168 a; + 984 a; 2 +380 a; 3 +35 a; 4 > ' 
r 59-148 ai + 79a; 2 + 10 a; 3 13-41 a;+23 a; 2 +5 a: 3 31-122 ai + 71 a; 2 + 20 a: 3 



L 329-1168 a; + 984 a; 2 + 380 a; 3 +35 a; 4 ' 329-1168 a; + 984 a; 2 +380 a; 3 +35 a; 4 ' 329-1168 a; + 984 a; 2 + 380 x 3 +35 a: 4 ' 

-18(-l + a:) 2 -i -i 

329-1168 a; + 984 a; 2 +380 a; 3 +35 a; 4 J J 



We compare the processor time conditioned by different algorithms for com- 
puting the Moore-Penrose inverse of constant matrices in the next table. Test 
matrices are taken from [27], and considered in the partial case a = 1. The 
test matrix name we state in the first column . Processor times required by the 
standard MATHEMATICA function PseudoInverse[] (see [26]) are allocated in the 
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second column of the table. Results corresponding to function Partitioning[] 
from [19] are placed in the third column. Fourth column is filled by the results 
generated by using the Leverrier-Faddeev algorithm from [9] . Results produced 
by applying MATHEMATICA implementation of the algorithm from [5] are placed 
in the next column, and the last two columns are arranged for the MATHEMATICA 
and DELPHI implementation of Algorithm 2.1. We use R = A in MATHEMAT- 
ICA functions ModGinvCholesky[) and DELPHI function ^41234 to compute the 
Moore-Penrose inverse. For matrix dimensions above 20 x 20 an application of 
the function ModGinvCholesky[] gives the information: "Result for Inverse of 
badly conditioned matrix << 1 >> may contain significant numerical errors"! 
These cases arc marked by the sign '*' in the table. Also, the sign '-' denotes a 
long processor time needed for the computation. 



Test 


Math. 


Math. 


Math. 


Math. 


Math. 


Delphi 


matrix 


Pseudolnverse 


Partitioning 


Lev.Faddeev 


Courrieu 


Alg. 2.1 


Alg. 2.1 


S5 


0.079 


0.016 


0.001 


0.001 


0.001 


0.062 


S10 


0.031 


0.031 


0.001 


0.015 


0.015 


0.062 


S25 




0.125 


0.062 


0.047 


0.109 * 


0.062 


S50 




1.187 


2.516 


0.375 


0.687 * 


0.940 


S100 




9.204 


44.375 


2.297 


5.781 * 


1.850 


F5 


0.125 


0.031 


0.001 


0.001 


0.001 


0.047 


F10 


1.094 


0.016 


0.001 


0.015 


0.015 


0.047 


F25 




0.047 


0.156 


0.110 


0.250 * 


0.062 


F50 




0.485 


2.672 


0.703 


2.328 * 


0.940 


F100 




2.812 


42.844 


5.782 


17.594 * 


1.850 


A5 


0.25 


0.006 


0.001 


0.001 


0.001 


0.047 


A10 


1.344 


0.015 


0.001 


0.015 


0.015 


0.062 


A25 




0.063 


0.171 


0.093 


0.265 * 


0.062 


A50 




0.484 


2.766 


0.766 


2.218 * 


0.940 


A100 




2.750 


43.781 


5.844 


16.954 * 


1.850 



Table 1. Processor time in Seconds for constant matrices 



4 Conclusion 

We introduce an algorithm for computing {1,2,3}, {1, 2, 4}-inverses, {2,3}, 
{2, 4}-inverses of prescribed rank as well as for computing the Moore-Penrose 
inverse for one-variable rational matrices. Our method uses the representa- 
tions of {i, j, fc}-inverses based on the matrix product involving the Moore- 
Penrose inverse and factors of the full-rank Cholesky factorization from [6] . On 
the other hand, a large number of representations and algorithms are avail- 
able for computing generalized inverses of rational and/or polynomial matrices 
[7, 8, 9, 10, 11, 17, 18, 20, 19, 22, 23, 2]. But, generalized inverses in these 
papers are computed using the Leverrier-Faddeev algorithm and the Grevile's 
algorithm. The algorithm proposed in this paper is an extension of the paper [5] 
to various classes of {i,j, fc}-inverses and to rational matrices. When the input 
matrix is constant, in a certain case R = A, we get an algorithm for computing 
the Moore- Pernose inverse, alternative with respect to the algorithm introduced 
in [5]. 
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Introduced algorithm is implemented in two different programming lan- 
guages: MATHEMATICA and DELPHI. The implementation in DELPHI is appro- 
priate only for constant matrices. In the constant matrix case we compare 
processor time required by these implementations of Algorithm 2.1 with respect 
to standard MATHEMATICA function Pseudoinverse, implementation of Grevile's 
partitioning method, implementation of Leverrier-Faddeev algorithm and the 
MATHEMATICA implementation of the algorithm from [5]. 

Column 2 is a confirmation of the statement that the method used in MATH- 
EMATICA function Pseudoinverse is time consuming for large matrices. The 
results from columns 3, 4, 5 and 6 in Table 1 again confirm known fact that 
MATHEMATICA (and other symbolic packages) is not applicable for large scale 
test problems. Our numerical experience shows that the algorithm introduced 
in [5] is superior with respect to the Grevile's partitioning algorithm for test 
matrices of smaller dimensions. But, the algorithm from [5] is inferior with re- 
spect to partitioning method in the case when test matrices of relatively great 
order from [27] are used. Leverrier-Faddeev algorithm produces the best results 
for test matrices of small dimensions and the worst results for test matrices of 
greater dimensions. 

Algorithm 2.1 produces inferior results with respect to algorithm from [5] 
for matrix dimensions greater than 20 x 20. The reason is clear. Algorithm 
from [5] computes the Moore-Penrose inverse using the Cholesky factorization 
of the matrix products A* A or AA* . On the other side, Algorithm 2.1 factor- 
izes the matrix products (A* A)* (A* A) or (AA*)(AA*)* , which produce bigger 
numbers causing badly conditioned matrices. But, our method for computing 
the Moore-Penrose inverse arises from a general algorithm, which is limited by 
the application of symmetric positive matrices (R*A)*(R*A) or (AT*)(AT*)* . 

5 APPENDIX 

For the sake of completeness we present the MATHEMATICA and DELPHI code for 
the implementation of Algorithm 2.1. 

5.1 Mathematica code 

In the following function we implement the Cholesky factorization. 

Cholesky [A0_ , n_] : =Module [{A=A0 , i , k,m , L ,U} , 
L=Table [0 , {n} , {n}] ; 
For [k=l,k<=n,k++, 

L[[k,k]]=Sqrt [A[[k,k]]-Sum[L[[k,m]] ~2 , {m, 1 ,k-l}] ] ; 
For [i=k+l , i<=n , i++ , 

L[[i,k]]=(A[[i,k]]-Sum[L[[i,m]]L[[k,m]] , {m, 1 ,k-l}] ) /L [ [k,k] ] ] ] ; 
U = Transpose [L] ; 
Return [L] ] 

In the auxiliary function Adop[a, j] we drop the last n — j columns from 
the matrix a. This function is used for the elimination of last zero rows in the 
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matrix Transpose[a\. 

Adop [a_List , j_] : =Module [{m, n}, 
{m,n}=Dimensions [a] ; 

Return [Transpose [Drop [Transpose [a] , - (n-j ) ] ] ] ; ] 

Function GinvCholesky [A] implements the algorithm from [5] 

GinvCholesky [A0_List] : =Module [{m , n , trans , A=A0 , L , M , Y} , 
{m,n}=Dimensions [AO] ; trans=False; 
If[m<n, trans=True; A=A0 .Transpose [AO] ; n=m, 

A=Transpose [AO] .AO] ; 
L=Cholesky [A,n] ; L=Simplif y [Adop [L , Matr ixRank [AO] ] ] ; 
M=Inverse [Transpose [L] .L] ; 

If [trans , Y=Transpose [AO] .L.M.M. Transpose [L] , 

Y=L .M.M. Transpose [L] .Transpose [AO]] ; 
Return [Simplify [Y] ] ] 

Function ModGinvCholesky [A, R] implements Algorithm 2.1. 

ModGinvCholesky [A_List ,R_List] : = 

Module [{m, ml.nl ,n,rr ,trans=False,L,M,Y,G,Gl}, 
{m,n}=Dimensions [A] ; {ml ,nl}=Dimensions [R] ; 
If [n==nl , trans=True ; 

G=Simplify [A. Transpose [R] . Transpose [A. Transpose [R] ] ] ; n=m, 

G=Simplify [Transpose [Transpose [R] .A] . Transpose [R] .A]] ; 
L=Cholesky[G,n] ;L=Adop [L,Min[MatrixRank[A] , MatrixRank [R] ] ] ; 
M=Inverse [Transpose [L] .L] ; 

If [trans ,Y=Transpose [R] . Transpose [A. Transpose [R] ] . L .M. M. Transpose [L] , 

Y=L. M.M. Transpose [L] .Transpose [Transpose [R] .A] .Transpose [R] ] ; 
Return [Simplify [Y] ] ] 

5.2 Delphi code 

We present the main part of DELPHI code for computing A{i, j, /c}-inverses of a 
given constant matrix A. Elementary functions used in computations are: func- 
tion TransMatQ which computes the transpose matrix, function MatMatRQ 
for the matrix multiplication, function MatrixRankQ for computing the matrix 
rank, the function which generates the matrix consisting of first i columns of 
a given matrix, called FifstlColumnsQ, and the function InversionMQ used 
for the usual matrix inversion. These functions are not restated here. 
Cholesky factorization is implemented in the following function. 

procedure TForml . Cholesky (AO : matrix ;var CO :matrix;n: integer) ; 

var i, j ,p,q: integer ;s : extended; si : real ; 

begin 

For i:=l to n do 

For j:=l to n do C0[i,j]:=0; 
For p:=l to n do 
begin 

s:=0; 

for q:=l to p-1 do s : =s+C0 [p,q] *C0 [p,q] 
sl:=A0 [p,p]-s; 

if sl<0. 00000000001 then sl:=0; 
C0[p,p] :=Sqrt(sl) ; 
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if CO[p,p]<>0 then 
begin 

for i:=p+l to n do 

begin 
s:=0; 

for j:=l to p-1 do s : =s+C0 [i , j] *C0 [p, j] ; 
C0[i,p] : = (A0[i,p]-s)/C0[p,p] ; 
end; 
end; 

end; 
end; 

Function A1234 implements Algorithm 2.1. 

procedure TForml . A1234(A1 ,R1 :matrix ;var L:matrix); 
var trans : boolean; minrank: integer ; 

Yl,Ll,L2,Gl,G2,G3,G4,G5,G6,G7,G8:matrix; 
begin 

trans :=false; 
if n=nn then 
begin 

trans : =true ; 

TransMat(Rl,mm,nn,Gl) ; MatMatR(Al ,G1 ,m,n,mm,G2) ; 
TransMat(G2,m,mm,G3) ; MatMatR(G2,G3,m,mm,m,Yl) ; 
n:=m; 
end 

else begin 

TransMat (Rl, mm, nn,Gl) ; MatMatR(Gl , Al ,nn,m,n,G2) ; 
TransMat(G2,nn,n,G3) ; MatMatR(G3,G2,n,nn,n,Yl) ; 
end; 

Cholesky(Yl,Ll,n) ; 
minrank :=MatrixRank (LI, n) ; 
firstIColumn(Ll , minrank, Yl) ; 

TransMat (Yl ,n, minrank, Gl) ; MatMatR(Gl , LI .minrank, n,n,G2) ; 
InverseM(G2,n,L2) ; 
if trans then 
begin 

TransMat(Rl,mm,nn,Gl) ; MatMatR(Al , Gl ,m,nn,m,G2) ; 
TransMat (G2,m,m,G3) ; MatMatR(Gl , G3,nn,mm,m,G4) ; 
MatMatR(G4,Ll,nn,m,n,G5) ; MatMatR(G5,L2,nn,n,n,G6) ; 
MatMatR(G6,L2,nn,n,n,G7) ; TransMat(Ll,n,n,G8) ; 
MatMatR(G7,G8,nn,n,n,L) ; 
WriteY(L,nn,n) ; 
end 

else begin 

MatMatR(Ll,L2,n,n,n,Gl) ; MatMatR(Gl,L2,n,n,n,G2) ; 
TransMat (LI, n,n,G3) ; MatMatR(G2,G3,n,n,n,G4) ; 

TransMat (Rl, mm, nn,G5) ; MatMatR(G5, Al ,nn,mm,n, G6) ; 

TransMat(G6,nn,n,G7) ; MatMatR(G7,G5,n,nn,mm,G8) ; 

MatMatR(G4,G8,n,n,mm,L) ; 
WriteY(L,n,mm) ; 
end; 

end; 
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